Real time evolution using the density matrix renormalization group 
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We describe an extension to the density matrix renormalization group method incorporating 
real time evolution into the algorithm. Its application to transport problems in systems out of 
equilibrium and frequency dependent correlation functions is discussed and illustrated in several 
examples. We simulate a scattering process in a spin chain which generates a spatially non-local 
entangled wavefunction. 

PACS numbers: PACS 
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The density matrix renormalization group (DMRG) Q 
is perhaps the most powerful method for simulating one 
dimensional quantum lattice systems. DMRG was origi- 
nally formulated as a ground state method. Later, it was 
generalized to give frequency dependent spectral func- 
tions [2|, |3( . The best spectral method, Jeckelmann's dy- 
namical DMRG Q, yields extremely accurate spectra. 
However, it is limited to only one momentum and one 
narrow frequency range at a time. Constructing an en- 
tire spectrum for a reasonable grid in momentum and 
frequency space can involve hundreds of runs. 

An alternative approach to dynamics with DMRG is 
via a real time simulation. Cazalilla and Marston intro- 
duced a real-time DMRG and used it to calculate the 
time evolution of one dimensional systems under an ap- 
plied bias[j|. In their approach, the DMRG algorithm 
is only used to calculate the ground state, and the time 
evolution in obtained by integrating the time-dependent 
Schrodinger equation in a fixed basis. Consequently, one 
expects it to lose accuracy when the wavefunction starts 
to differ significantly from the ground state. In the sys- 
tems studied by Cazalilla and Marston, the time evolu- 
tion could be carried out for a reasonable length of time 
before this happened. Luo, Xiang, and Wang 6] showed 
how to construct a basis which applies to a time-evolving 
wavefunction over a whole range of times simultaneously. 
This approach was shown to be more accurate than that 
of Cazalilla and Marston, but it is not very efficient — the 
basis must be quite large to apply to a long interval of 
time, and the whole time evolution is performed at every 
DMRG step. 

Recently, Vidal developed a novel time-dependent sim- 
ulation method for near-neighbor one dimensional sys- 
tems which overlaps strongly with DMRG 7J. The 
crucial new idea of the method is to use the Suzuki- 
Trotter decomposition for a small time evolution operator 
exp(— itH). The second order Suzuki- Trotter break-up 
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terms exp(—irHj) (coupling sites j and j + 1) within Ha 
or Hb commute. Writing the wave function in matrix 
product form (which underlies the DMRG block form 8J), 
Vidal showed that one can apply each link-term directly 
to the wavefunction, exactly and efficiently. After each 
such application, a Schmidt decomposition, equivalent to 
diagonalizing the DMRG density matrix, is performed to 
return the wavefunction to the matrix product form. One 
applies all the Ha terms, and then all the Hb terms, etc. 

Although this method seems very efficient, a number 
of aspects are novel for DMRG users, stemming from the 
fact that one does not ordinarily deal with the matrix 
product representation directly. Implementing this idea 
into a DMRG algorithm may be time consuming and may 
require a very substantial rewriting of one's program. 

In this paper we take the key idea of the Suzuki- Trotter 
decomposition, but we modify it and apply it in a more 
natural way within the context of DMRG. The result is 
a surprisingly simple yet very powerful modification of 
the algorithm for real-time dynamics which we believe 
can be incorporated into a typical program in only a day 
or two of programming. We illustrate the approach with 
real-time simulations which set a new paradigm for the 
size and accuracy obtainable. 

The standard DMRG representation of the wavefunc- 
tion at a particular step j during a finite system sweep 
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Here we have a left block containing many sites (with 
states I), two center sites (with states ay, aj+i), and a 
right block (states r). The states I and r are formed as 
eigenvectors of a density matrix, and represent a highly 
truncated but extremely efficient basis for representing 
the ground state, plus any other targeted states which 
have been included in the density matrix. Now suppose 
we have an arbitrary operator A acting only on sites j 
and j + 1. This operator can be applied to \ip) exactly, 
and reexpressed in terms of the same optimal bases, with 



where Ha contains the terms of the Hamiltonian for the 
even links, and Hb for the odd. The individual link- 
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If A included terms for other sites, we could not write 
this simple exact relation; new bases would need to be 
adapted to describe both \tp) and A\tp), requiring perhaps 
several finite-system sweeps through the lattice. 

This implies that we can apply the link time evolution 
operator exp(— irHj) exactly on DMRG step j, but at 
any other DMRG step the application is approximate. 
Accordingly, we adapt the Suzuki- Trotter decomposition 
to match the DMRG finite-system sweeps, so that each 
term can be applied exactly. We decompose the time 
propagator as 

e -irH _ e -irH 1 /2 e -irH 2 /2 _ _ _ e -,rH 2 /2 ( ,-,rH 1 /2^ /j\ 

This decomposition is good to the same order (errors 
of order r 3 ) as the usual odd/even link decomposition. 
When applied 1/r times to evolve one unit of time, the er- 
rors are r 2 . The main idea is then to apply exp(— irHi/2) 
at DMRG step 1, then exp(— irif 2 /2) at step 2, etc., 
forming the usual left-to-right sweep, then reverse, apply- 
ing all the reverse order terms in the right-to-left sweep. 

This procedure requires one to use the step-to- 
step wavefunction transformation first developed to 
provide a good guess for the Lanczos or Davidson 
diagonalization 9] , where it can improve run times by an 
order of magnitude. It transforms the wavefunction from 
the basis of step j ± 1 to that of step j. Assuming this 
transformation is present in a ground state DMRG pro- 
gram, the real-time algorithm introduces only a very mi- 
nor modification: at step j, instead of using the Davidson 
method, one evolves the transformed wavefunction by ap- 
plying exp(— irHj/2). Before the time evolution starts, 
we typically use ordinary DMRG to find the ground state. 
Next, we either (1) change the Hamiltonian, or (2) apply 
an operator to the ground state to study a new wave- 
function which is a combination of excited states. 

As a first test case, of type (1), we study the models of 
Eqs.(2) and (5) in Ref. |5j corresponding to a quantum 
dot connected to two non-interacting leads, and a junc- 
tion between two Luttinger liquids, respectively, driven 
out of equilibrium by a voltage bias. In these cases at 
t = a bias in the chemical potential is turned on as 
a smoothed step function, making the new Hamiltonian 
time-dependent. At each time step, the expectation value 
of the current operator (defined by Eq.(4) of Ref.[j|) is 
calculated. In FigQ]we show the results for a chain of 
length L = 64 and the same set of parameters used in 
Ref.[||, keeping only m = 128 states and using a time 
step r = 0.2. It should be compared to Figs. 1 and 2 
in Ref. || and Figs. 1 and 2 in Ref.0- Our results 
exceed the accuracy obtained by the previous methods, 
with only a fraction of the states. For the non-interacting 
problem, the agreement with the exact results is excel- 
lent up to times t ~ 70. The main reason for obtaining 
higher accuracy for fixed m compared to Ref. is that at 
any step we only need to target one state at one instant 
of time, not a whole range of times. 
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FIG. 1: (a) Tunneling current through a non-interacting 
quantum dot and a junction as defined in Eqs. (2) and (5) 
in Ref.0, respectively. The full lines correspond to exact re- 
sults, (b) Tunneling current through an interacting junction, 
with V = 0.5 and V = 1.1. All the DMRG results where 
obtained with M = 128 and a time step r = 0.2. 



In method (2), after finding the ground state \(f>), we 
apply an operator A at t = 0, to obtain \ip(t — 0)}, 
and evolve in time. We first use this approach to calcu- 
late time-dependent correlation functions. In this case, 
we time evolve both \4>(t)) and \ip(t)), including both as 
target states for the DMRG density matrix. Although 
the time dependence exp(-iEot) of \4>(t)} is known {Eq 
is the ground state energy), by evolving it we keep its 
representation in the current basis. In addition, we ex- 
pect a significant cancelation in the errors due to the 
Suzuki- Trotter decomposition in constructing the corre- 
lation functions. A typical correlation function is calcu- 
lated as 

< <p\B(t)A(0)\<fi >=< m\B\m >, (5) 

We use a complete half-sweep to apply A to \4>). In par- 
ticular, if A is a sum of terms Aj over a number of sites, 
then we apply an Aj only when j is one of the two cen- 
tral, untruncated sites. Thus the basis is automatically 
suitable for Aj\<j>). During this buildup of A at step j we 
target both the ground state \<f>) and J2j'=i AH^)- 

As an example we consider the spin-1 Heisenberg 
chain, with Hamiltonian 

H = J2SjS j+1 , (6) 

3 

where we have set the exchange coupling J to unity. This 
system has a gap (the "Haldane gap") of A# = 0.4105 to 
the lowest excitations, which are spin-1 magnons at mo- 
mentum 7r, and a finite correlation length of £ = 6.03. |l0j 
The single-magnon dispersion relation has been calcu- 
lated with excellent accuracy 3]. However, determina- 
tion of the full magnon line is quite tedious with existing 
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FIG. 2: Time evolution of the local magnetization (S z (x)) of 
a 200 site spin-1 Heisenberg chain after S' + (100) is applied. 



DMRG methods. Here we demonstrate how to calculate 
the entire magnon spectrum with only one time depen- 
dent DMRG run. 

We take A = S + (j) for a single site j in the center 
of a long chain. This operator constructs a localized 
wavepacket consisting of all wavevectors. This packet 
spreads out as time progresses, with different components 
moving at different speeds. The speed of a component is 
its group velocity, determined as the slope of the disper- 
sion curve at k. In Fig[21we show the local magnetiza- 
tion (ip(t)\S z \tjj(t)) for a chain of length L = 200, with 
timestep r = 0.1. At t = 0, the wavepacket has a fi- 
nite extent, with size given by the spin-spin correlation 
length £. At later times, the different speeds of the dif- 
ferent components give the irregular oscillations in the 
center of the packet. We kept m — 150 states per block, 
giving a truncation error of about 6 x 10~ 6 . 

^From this type of simulation we can construct the 
Green's function 



G(x,t) = -i(<t>\T[S^(t)S+(0)}\<t>) 



(7) 



as G(x,t) — — i((/)(\t\)\S~\ijj(\t\)). Here x is measured 
from the site j where S + is applied. We make one mea- 
surement of G(x,t) for each left-to-right DMRG step, 
namely for step x. For efficiency we measure as we evolve 
in time, rather than, say, devoting every other sweep to 
measurements without time evolving. This may worsen 
the Suzuki- Trotter error somewhat, but we have found 
the results quite satisfactory. Since G(x, t) is even in x 
and t, the Fourier transform is 
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cos kxG(x, t) (8) 



The spectral function of interest is — l/7rImG(fc, u>). We 
inevitably have some cutoff in time T for the available 
data. The maximum useable value of T depends on 
the length of the chain: when the leading edge of the 




FIG. 3: The single magnon line of the spin-1 Heisenberg an- 
tiferromagnetic chain. The entire spectrum is obtained from 
one DMRG run, by Fourier transforming the time and posi- 
tion dependent correlation function {Sf (t)Sg (0)}. The broad 
solid curve shows the location of the maximum in the spectra 
for a particular q, in units of the Haldane gap, 0.41050(2), for 
a system of L = 600 sites, using a time step r = 0.02, running 
for T = 27.3, and keeping m = 200 states. For comparison, 
results from two other runs are shown: L — 400, r = 0.1, 
T = 60, and m = 150 (dashed curve); and. L = 400, r = 0.4, 
T — 72, and m — 200 (dotted curve). The solid curve peaked 
at q — it, shown only for the first run, is the weight Ao in this 
quasiparticle peak, i.e. S(u) « Aq5(u} — luq). 



wavefront hits the ends of the chain, the data no longer 
describes an infinite chain. On the other hand, before 
that point the data does describe an infinite chain with 
boundary effects dying off exponentially from the edges. 
This allows us to precisely specify the momentum k, for 
times t < T. To perform the time integral we multi- 
ply G(x, t) by a windowing function W(t) which goes 
smoothly to zero as t — ► T. We have chosen a Gaussian, 
exp(— 4(t/T) 2 ), which has the advantage of having a non- 
negative Fourier transform, yielding a nonnegative spec- 
tral function (except for possible terms of size exp(— 4)). 
Note that if the true spectral function has an isolated 
delta function peak, the windowed spectrum will have a 
Gaussian peak centered precisely at the same frequency. 
Thus it is possible to locate the single magnon line with 
an accuracy much better than l/T. If a continuum is 
also present nearby, the peak is less well determined. In 
the case of the S — 1 chain, for k near n the peak is iso- 
lated, but at some point near k — 0.37T the peak enters 
the two magnon continuum and develops a finite width. 
Note that from our single simulation we determine the 
spectral function for a continuum of values of k and u>. 

In FigEl we summarize the results for the single 
magnon peak, determined automatically as the maxi- 
mum of the spectrum. To gauge the errors we present 
results for several runs with various parameters. The re- 
sults are very close; the most signifcant errors are due 
to a finite r, with the r = 0.4 run showing some non- 
negligible errors. The results agree very well with ac- 
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FIG. 4: A gaussian magnon wavepacket with momentum 
k = 0.8-7T, S z — 1, and halfwidth 16 scattering off the left end 
of a 200 site spin-1 chain. The chain end states initially have 
S z — —1/2. In all cases we measure (S z (x)} at time t, and 
r = 0.2. In (a) we show t = and t = 20. In (b), we show the 
left-most 50 sites for a number of times. In (c), we show the 
whole chain for t = and t — 110. After the scattering, the 
system is in a non-local superposition of spin-flip and non- 
spin-flip states. 



either with or without a spin flip occuring. If the spin 
flip occurs, the end spin changes to S z = 1/2 and the 
wavepacket to S z = 0. We see from the figure that af- 
ter the scattering, the end spin seems to have taken on 
an intermediate value of S z , in particular (S z ) « —0.11. 
Meanwhile, the wavepacket seems to have a total spin of 
(S z ) s» 0.61. Angular momenta of each are still quan- 
tized. The intermediate values occur because we are ob- 
serving a "macroscopic" quantum superposition of the 
state with and without the spin flip. Specifically, the 
scattering is described by 



|,1> — a|— |,l>+b||,0) 



(9) 



where a 2 w 0.61, b 2 rj 0.39. Note that the scattered 
S z = magnon does not show up when we measure 
S z . A local description of the wavefunction, as (a\ — ^} + 
x (t 1 1) +5|0)), is not possible; it does not conserve 
total S z . Our measurement of S z , {ip(t)\S*\ip(t)), does 
not affect the state, but if one performed a real exper- 
iment and measured the spin of, say, the magnon after 
scattering one would obtain either S z = 1 or S z = 0. 
We find it quite remarkable that we can simulate a pro- 
cess in which such a non-local superposition develops! 
Our results raise the possibility of using such spin chains 
for experimental studies of quantum measurement and 
quantum computation. 

We acknowledge the support of the NSF under grants 
DMR-0311843. 



curate frequ ency -based DMRG results 3] and quantum 
Monte CarloHJ. 

With real-time dynamics, we can simulate processes 
which would be very difficult to understand via fre- 
quency dynamics. As an example, we consider a magnon 
wavepacket scattering off the end of a spin-1 chain, shown 
in Fig0J The magnon is a triplet, with S z = 1, and trav- 
els to the left with a speed of about 2.0. The open ends 
of a spin-1 chain have spin-1/2 degrees of freedom, which 
have received considerable attention [ToL Il2| . One can 
view this end state as a spinon bound to the end. An an- 
tiferromagnetic oscillation accompanies this state, decay- 
ing exponentially with the correlation length away from 
the edge. We choose the ground state with total spin 
S z = —1, making the end states each have S z = —1/2. 
When the wavepacket hits the left end, it can scatter 



[1] S.R. White, Phys. Rev. Lett. 69, 2863 (1992), Phys. Rev. 

B 48, 10345 (1993). 
[2] K.A. Hallberg, Phys. Rev. B 52, 9827 (1995). 
[3] T.D. Kuhner and S R. White, Phys. Rev. B 60, 335 

(1999). 

[4] E. Jeckelmann, Phys. Rev. B 66, 045114 (2002). 

[5] M. A. Cazalilla and J. B. Marston, Phys. Rev. Lett. 88, 

256403 (2002); see also 91, 049702 (2003). 
[6] H. G. Luo, T. Xiang, and X. Q. Wang, Phys. Rev. Lett. 

91, 049701 (2003). 
[71 G. Vidal, P hys. Rev. Lett. 91, 147902(2003); and 

|quant-ph /0310089 
[8] S. Ostlund and S. Rommer, Phys. Rev. Lett. 75, 3537 

(1995). 

[9] S.R. White, Phys. Rev. Lett. 77, 3633 (1996). 
[10] S.R. White and D.A. Huse, Phys. Rev. B 48, 3844 (1993). 
[11] M. Takahashi, Phys. Rev. Lett. 62, 2313 (1989); Phys. 

Rev. B 48, 311 (1993). 
[12] I. Affleck, T. Kennedy, E. H. Lieb and H. Tasaki, Phys. 

Rev. Lett. 59, 799 (1987). 



